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While universal quantum computers ideally 
solve problems such as factoring integers ex- 
ponentially more efficiently than classical ma- 
chines, the formidable challenges in building 
such devices motivate the demonstration of 
simpler, problem-specific algorithms that still 
promise a quantum speedup. We construct a 
quantum boson sampling machine (QBSM) to 
sample the output distribution resulting from 
the nonclassical interference of photons in an 
integrated photonic circuit, a problem thought 
to be exponentially hard to solve classically. 
Unlike universal quantum computation, boson 
sampling merely requires indistinguishable pho- 
tons, linear state evolution, and detectors. We 
benchmark our QBSM with three and four pho- 
tons and analyze sources of sampling inaccu- 
racy. Our studies pave the way to larger devices 
that could offer the first definitive quantum- 
enhanced computation. 

Universal quantum computers require physical 
systems that are well-isolated from the decohering 
effects of their environment, while at the same time 
allowing precise manipulation during computation. 
They also require qubit-specific state initialization, 
measurement, and the generation of quantum cor- 
relations across the system (1-5). Although there 
has been substantial progress in proof-of -principle 
demonstrations of quantum computation (6-10), si- 
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Figure 1 : Model of quantum boson sampling. Given a spec- 
ified initial number state jT) = \T\...Tm) and linear trans- 
formation A, a quantum boson sampling machine efficiently 
samples from the distribution P(S|T) of possible outcomes 
|S> = \Si...S M ). 



multaneously meeting these demands has proven 
difficult. This motivates the search for schemes that 
can demonstrate quantum-enhanced computation 
under more favorable experimental conditions. The 
power of one qubit (11), permutational (12) and in- 
stantaneous (13) quantum computation, are exam- 
ples that have allowed investigation of the space be- 
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Figure 2: Schematic and characterization of the photonic 
circuit, (a) The silica-on-silicon waveguide circuits consist 
of M = 6 accessible spatial modes (labeled 1—6). For the 
three-photon experiment, we launch single photons into in- 
puts i = 2, 3 and 5 from two parametric down-conversion 
sources and postselect outcomes in which three detections are 
registered amongst the output modes j. For the four-photon 
experiment, which is implemented on a different chip of iden- 
tical topology, we inject a double photon pair from a single 
source into the modes i = 1, 3 and postselect on four detec- 
tion events, (b)-(c) Measured elements of the linear transfor- 
mation Aij = Ttje" l ' ij linking the input mode i to the output 
mode j of our three-photon apparatus. The circuit topology 
dictates that several m are zero, and our phase-insensitive in- 
put states and detection methods imply only six non-zero 4>ij ■ 
Since only relative values are needed due to post-selection, we 
rescale each row of r so that its maximum value is unity. 

tween classical and universal quantum machines. 

Boson sampling has recently been proposed as 
a specific quantum computation that is more ef- 
ficient than its classical counterpart but only re- 
quires identical bosons, linear evolution, and mea- 
surement (14). The distribution of bosons that 
have passed through a linear system represented 
by the transformation A is thought to be exponen- 
tially hard to sample from classically (14). The 
probability amplitude of obtaining a certain out- 
put is directly proportional to the permanent of 
a corresponding submatrix of A. The permanent 
expresses the wavefunction of identical bosons, 
which are symmetric under exchange (15, 16); 
in contrast, the Slater determinant expresses the 



wavefunction of identical fermions, which are an- 
tisymmetric under exchange. While determinants 
can be evaluated efficiently, permanents have long 
been believed to be hard to compute (17); the best 
known algorithm scales exponentially with the size 
of the matrix (18). One can envision a race be- 
tween a classical and a quantum machine to sample 
the boson distribution given an input state and A. 
The classical machine would evaluate at least part 
of the probability distribution, which requires the 
computation of matrix permanents. The QBSM in- 
stead creates indistinguishable bosons, directly im- 
plements A, and records the outputs. For a suffi- 
ciently large problem size, the quantum machine is 
expected to win (14). 

Photonics is a natural platform to implement 
boson sampling since sources of indistinguishable 
photons are well-developed (19), and integrated 
optics offers a scalable route to large linear sys- 
tems (20). In addition, photons have extremely 
low decoherence rates. Importantly, boson sam- 
pling requires neither nonlinearities nor on-demand 
entanglement, unlike photonic approaches to uni- 
versal quantum computation (21). This clears the 
way for experimental boson sampling with existing 
photonic technology, building on the extensively 
studied two-photon Hong-Ou-Mandel (HOM) in- 
terference effect (22). 

A QBSM (Fig. 1) samples the output distribution 
of a multi -particle bosonic quantum state |\P ut)j 
prepared from a specified initial state |T) and lin- 
ear transformation A. A trial begins with the input 
state |T) = \Tl.Tm) oc rL=i(«!) Tl |0}, which 
describes N = ^ particles distributed in M 

input modes in the occupation-number representa- 
tion. The output state l^out) is generated according 
to the linear map from output to input mode cre- 
ation operators a\ = 2~2j=i Ajj&j, where A is an 
MxM matrix. Finally, the particles in each of the 
M output modes are counted. The probability of a 
particular measurement outcome |S) = \S\...Sm) 
is given by 

|Per(A( S ' T ) N )| 2 

p(s|t) = i<si* out >i 2 = — l i ;\ 

llj=i °i- lli=i L i- 

(1) 

where the NxN submatrix A( S ' T - ) is obtained by 
keeping Sj (T{) copies of the j th column (i th row) 
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Figure 3: Boson sampling results. The measured relative frequencies P cxp of outcomes in which the photons are detected in 
distinct modes are shown in red for (a) three- and (b) four-photon experiments. Each data set is collected over 160 hours, and 
statistical variations in counts are shown by the red shaded bars. The theoretical distributions P th (blue) are obtained from the 
permanent of a submatrix constructed from the full transformation A, as depicted in the inset. The blue error bars arise from 
uncertainties in the characterization of A. 



of A (23). 

The aim of boson sampling is to sample from the 
distribution given by P(S|T). A QBSM achieves 
this efficiently by directly sampling KSj^out)) 2 - 
The classical strategy, on the other hand, requires 
an estimation of the right-hand side of Eqn. (1), for 
which the calculation time scales exponentially in 
N(18). 

Our QBSM consists of sources of indistinguish- 
able single photons, a multiport linear optical cir- 
cuit, and single-photon counting detectors. Two 
parametric down-conversion (PDC) pair sources 
(24) are used to inject up to four photons into a 
silica-on-silicon integrated photonic circuit, fabri- 
cated by UV writing (25). The circuit is shown in 



Fig. 2a and consists of M=6 input and output spa- 
tial modes coupled by a network of ten beam split- 
ters (20). The output state is measured with single- 
photon avalanche photodiodes on each mode. We 
only consider outcomes in which the number of de- 
tections equals the intended number of input pho- 
tons (23). 



Our central result of three- and four-boson sam- 
pling is shown in Fig. 3. In the first case, we repeat- 
edly inject three photons in the input state |T) = 
|011010), monitor all outputs, and collect all three- 
fold coincident events. In the four-photon experi- 
ment, we use the input | T) = 1 202000) and record 
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Figure 4: Sampling accuracy. We consider several boson outcome distributions: the experimental samples p cx P ; the ideal 
predictions of Eqn. (1) P th , and the predictions of the full model p mod . (a) The L\ distance d between P th and a Monte Carlo 
simulation of an ideal machine that samples P th a finite number of times for our three and (b) four photon cases. The inset 
histograms show the variation in d expected for a sample size corresponding to the 1421 and 405 counts collected in our three- 
and four-photon experiments. The distance d(P oxp ,P th ) (red) between our experimental sample and the ideal distribution 
suggests an underlying systematic inaccuracy. A full model which includes multiphoton emission and photon distinguishability 
is validated by the distance d(P cxp , p mod ) (green) which is consistent with statistical variation, (c) The predicted variation 
in d(P th , p mod ) i s shown as a function of A for the photon distinguishabilites V indicated for the three- and (d) four-photon 
cases. Our experimental conditions are marked (red dot). 



all four-fold events. For each experiment, the 
measured relative frequencies P$ xp for every al- 
lowed outcome |S) are shown along with their ob- 
served statistical variation. The corresponding the- 
oretical Pg h , calculated using the right-hand side of 
Eqn. (1), are shown along with their uncertainties 
arising from the characterization of A, described 
below. 

We reconstruct A with a series of one- and two- 
photon transmission measurements to determine its 
complex- valued elements Ay = r^e^ (26). The 
characterization results for the circuit used in the 
three-photon experiment are shown in Fig. 2b-c. 
To obtain the magnitude T y -, single photons are 
injected in mode i. The probability of a subse- 
quent detection in mode j is given by P\(j,i) = 
|Ajj| 2 = r? . The phases faj are determined 

'Eqn. (1) is expected to hold for any jS) and |T), how- 
ever the classical hardness of sampling P(SjT) is maximized 
when Sj,Ti G {0,1} for a given N and A. 



from two-photon quantum interference measure- 
ments. The probability that a photon is detected 
in each of modes j\ and j% when they are injected 
in modes i\ and ii is given by P<2.{j\, ji,i\,i%) = 
Ai 2 j 2 + Ai 2 j 1 A.i 1 j 2 | 2 . This expression is used 
to relate the relevant phases 0jj given the previ- 
ously determined magnitudes T{j. We use a least- 
squares fit to determine the phases 4>ij from a com- 
plete set of two-photon measurements (23). 

To analyze the performance of our QBSM we 
compare our results to an ideal machine. We 
quantify the match of two sets of relative fre- 
quencies P^ 1 ) and P^ 2 ) by calculating the L\ dis- 
tance dW(P«,p( 2 ))=± Zs where 
N denotes the number of photons in a sample (27). 
Identical and maximally dissimilar distributions 
correspond to d=0 and d=l, respectively. For 
our experiments we calculate S N \P cxp , P th ) to 
give d< 3 )=0.094 ± 0.014 and =0.097 ± 0.004. 
Even in an ideal QBSM with perfect state prepa- 
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ration and detection, the statistical variations result 
in nonzero d. If we substitute for our experimental 
data a Monte Carlo sampling of P th with sample 
size equivalent to our experiments, we instead cal- 
culate d( 3 )=0.043±0.012 and (i( 4 )=0.059±0.022. 
This suggests there are significant contributions to 
d(P exp , P th ) beyond statistical deviation. 

Due to experimental limitations, our QBSM oc- 
casionally samples distributions other than P th . 
The dominant sources of this sampling inaccuracy 
are multi-photon emission and partial distinguisha- 
bility of our photon sources. In practice, all single- 
photon sources produce multiple photons with a fi- 
nite probability (19). For our PDC sources, the out- 
put state is approximately |00) + A|ll) + A 2 |22), 
with A^l. Both single -photon and undesired mul- 
tiphoton terms increase with A. In our three-photon 
experiments, for example, multiphoton emission 
from the two PDC sources can lead to input states 
|T) = |021010) or |012020), which contribute 
to three-fold coincident events if photons are lost 
or emerge in the same output mode. In addition, 
partial distinguishability of the photons contami- 
nates the distribution sampled by the QBSM by 
mixing in one- and two-photon interference effects 
(28). To characterize the photon indistinguisha- 
bility we performed two-photon Hong-Ou-Mandel 
experiments, in which interference visibilities were 
90% of the expected values on average. 

We form a new distribution p mod that ac- 
counts for the effects of multiphoton emission 
and photon distinguishability (23). The distance 
d(P ex P, pmod^ s h own by the green point in the in- 
sets of Fig. 4a-b, is found to be consistent with the 
statistical variation due to a finite sample size, for 
both the three- and four-photon experiments. This 
suggests we have correctly identified and modeled 
the significant sources of inaccuracy. Our validated 
full model may be used to predict the performance 
of realistic QBSMs and guide the design of future 
machines with larger N. To investigate how the 
performance of our QBSM depends on A and pho- 
ton distinguishability, we calculate <i(P mod , P th ) 
for a range of operating parameters (Fig. 4c-d). In 
terms of A, a clear tradeoff is presented between 
data rate and inaccuracy due to multiphoton emis- 
sion, which is an intrinsic consequence of using 
PDC sources. Improvement in photon indistin- 



guishability increases the fidelity to the ideal ma- 
chine, and additionally is thought to enhance the 
computational power of a QBSM (28). 

Our results demonstrate that boson sampling 
provides a model of quantum-enhanced computa- 
tion that is experimentally feasible with existing 
photonic technology. Future generations of QB- 
SMs will benefit from ongoing advances in in- 
tegrated photonics such as reduced transmission 
loss, efficient number-resolving detectors (29), and 
multiplexed (30, 31) or single-emitter (19) photon 
sources. Such improvements suggest a promising 
route to large-scale photonic QBSMs, which will 
provide clear evidence for the computational power 
of quantum mechanics. 
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1 Supplementary Information 

In this section we first outline how the boson distri- 
bution is given by a set of matrix permanents. We 
then show how the boson distribution can be accu- 
rately predicted by the permanents associated with 
a non-unitary matrix describing a lossy channel, 
A, if one post-selects on trials where no photons 
are lost. We describe how one- and two-photon 
statistics can accurately characterize A. Finally, 
the principal sources of error in this experiment, 
namely the photon distinguishability and higher or- 
der terms from our PDC photon sources, are dis- 
cussed. 



1.1 The boson distribution is given by a set 
of matrix permanents 

We assume N photons are injected into an M mode 
circuit. We will consider the simplified case where 
the input (and output) states contain no more than 
one photon per mode, though the general case is 
treated elsewhere (32). Without loss of generality, 
let modes 1 to N contain an input photon, while 
modes N + 1 to M have vacuum inputs. The input 
state can then be described by 



N 



i*in>=iT>=n«!io> 



(2) 



8=1 



The unitary transformation allows one to evolve the 
operators according to 



M 



(3) 



i=i 



where a\ and bj are creation operators on the i-th 
input and j-th output mode respectively. We then 
obtain the output state 



N 



M 



n em i°) 



(4) 



=i \j=i 



To find the boson distribution, we project our 
output onto a state |S) which, in the number state 
basis, we describe by an N element vector S, where 



Sj gives the mode of the j-th photon. The proba- 
bility of measuring this state is then 



P5=|(Sh/Wt)f 



(5) 



N 



N 



M 



3=1 \k=l 



The term in square brackets can be expanded and 
includes M N terms, as one is selecting N photons 
from M modes where repetitions are allowed (> 
1 boson in a mode). One can rewrite this term in 
square brackets to give 



Ps 



N 



<oin s * 



N 



j=i \k=i 



|0) 



(6) 

where V is the set of M N permutations of N pho- 
tons amongst M modes, repetitions allowed. The 
tilde notation will be used throughout this paper for 
a set of permutations. Then, indicates the spa- 
tial mode of the A;-th photon in the j-th permuta- 
tion. As an example, consider the case with M = 3 
modes and N = 2 input photons, then 



V 1 = [l,l] V 2 
V 4 = [2,l] V 5 
V 7 = [3, 1] V s 



[1,2] V 3 
[2, 2] V 6 
[3,2] V 9 



[1,3] 
[2,3] 
[3,3] 



(V) 



Let us denote all AH permutations of S by S where 
S J k indicates the spatial mode of the k-th photon in 
the j-th permutation. For example, if we project 
onto the state \S) = |011) then S = [2, 3] and 



S 1 = [2, 3] S 2 = [3, 2] 



(8) 



It is clear we only retain terms from the summation 
in Eqn. (6) where G S, otherwise at least one 
annihilation operator will act on vacuum and give 
Ps = 0. This then leaves us with 

2 



Ps 



N\ N 

£IRs; 

j=l k=l 



(9) 



The formula for the permanent of an n x n matrix 

A = (a,ij) is 



Per(^) 



n\ n 
i=l j=l 



(10) 
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where <t* gives the j-th element of the i-th permu- 
tation of the numbers 1, 2, n. The term inside 
the modulus of Eqn. (9) has the same form as the 
matrix permanent in Eqn. (10). Our original uni- 
tary, U, can be described by an M x M matrix. 
However, it is obvious from Eqn. (9) that, in gen- 
eral, we take the permanent of a subsection of U. 
Specifically, we only keep rows 1 — >• N, those rows 
corresponding to modes with input photons. In ad- 
dition, we only keep columns corresponding to the 
elements of S. Let us call this modified subsection 
of our original unitary U^ S ' T \ Then, using the def- 
inition of the permanent we can rewrite Eqn. (9) 



P(S|T) = Per(?7 (S ' T) ) 



(11) 



If one allows the possibility of more than one pho- 
ton per input and output mode, then a similar anal- 
ysis yields Eqn. (1) (32). 

1.2 Effects of loss 

In any experimental implementation of boson sam- 
pling there will be losses. These losses, regardless 
of where they occur, can be modeled as beam split- 
ters that link accessible to inaccessible modes (33). 
When these losses are considered, it is important to 
ask whether the boson distribution is still given by 
a set of matrix permanents and, if so, what linear 
transformation does that matrix describe. 

Let us adopt the convention that modes 1 to M 
are accessible modes while loss modes are given 
the labels M + 1 to L. There is an L x L unitary 
operation describing the evolution of our pure input 
state, though we must trace over these loss modes 
at the output, yielding a mixed state over the acces- 
sible modes. Let us assume that photons are input 
into modes 1 to TV where N < M and we post- 
select on cases where N photons are detected, by 
definition, in the accessible modes 1 to M. Then 
Eqn. (6) becomes 



Ps 



N 




(12) 

but Si < M as we can only project on accessi- 
ble modes. Therefore, even though U is an L x L 
matrix and the elements of V range from 1 — >• L, 



when we project onto S (all the permutations of S) 
we are left with 



Ps 



N\ N 



3=1 k=l 



k,Si 



Per(u( s < T )) = Per(A( s ' T ) 



(13) 



where [7^ S ' T - ) is again a modified version of the 
original unitary but only keeping rows 1 — > N and 
columns in S[ < M. Since these elements always 
describe the accessible modes, then we can equiv- 
alently work in terms of A, where A*j = Uij but 
i,j < M. 

In summary, when post-selecting on no bosons 
being lost, one can work in terms of A, a non- 
unitary linear transformation that is simply the sub- 
section of U over the accessible modes. 

Equivalently, we can describe our system as a 
noisy (lossy) quantum channel in the operator sum 
representation (34). This formalism will be useful 
later when we discuss sources of error from higher 
order PDC terms. In this picture, the accessible 
modes are in the system Q, while all inaccessible 
loss modes form the environment system E. We 
can describe the transformation induced by our cir- 
cuit over the full space Q ® E by a unitary oper- 
ation U. Let p and a be the inputs to Q and E 
respectively, then the output in Q after a projective 
measurement P m and tracing over the environment 
is described by 



Pout 



\XE(PmU(p®a)U^P m ) (14) 



Let the basis for E be described by je^) and the 
initial state of the environment be a = Qj\j) {j \ . 
then we can express Eqn. (14) as 



Pout = Yj E 3kP E jk 



(15) 



where Ejf. = ^/qj(ek\P m U\j) are the Kraus op- 
erators. We do not directly characterize U, as it 
extends over the environment which is inaccessible 
to the experimenter. However, with photons we can 
assume a = |0)(0|, and all of our boson sampling 
results post-select on the case where no photons are 
lost to the environment. Therefore, the summation 
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in Eqn. (15) reduces to only one term with postse- 
lection 



Pout — EqqpE { 



do 



(16) 



Experimental boson sampling efforts will sample 
a non-unitary transformation that is equivalent, 
when post-selecting on no bosons being lost, to the 
Eqq Kraus operator. 



1.3 Characterizing the unitary over the 
accessible modes 

We want to find Ay, the elements of the non- 
unitary linear transformation induced by our lossy 
circuit, to predict the post-selected boson distri- 
bution. Every element is a complex number that 
can be written Ay = rye*^. We can find r^jj 
by coupling light into mode i\ and monitoring the 
probability of output in j\ yielding 



MIjo> 



M 



(OI&^A^IO) 

3=1 



\A 



'nji i 



(17) 



It is sufficient to measure the ratio of values be- 
tween output ports, {Ti 1 j 1 /ri 1 j 2 ) 2 , as this de- 
scribes the transformation induced by our circuit 
up to a constant factor for each input, Xj. If A is 
the true transformation, then this procedure gives 
us A A where A is a diagonal matrix with entries 
Xa = Xi. However, as Per (A A) = £j)Per(A) 
then Eqn. (1) shows, since we always launch the 
same input state, every portion of our boson distri- 
bution is multiplied by a constant factor that can- 
cels in the normalization. These data are collected 
by periodically pausing boson sampling data col- 
lection, blocking two of the three photon inputs, 
and measuring the relative power in each output 
mode. We thus obtain accurate values for r as well 
as a variance that is important for calculating the 
error bars in Fig. 3. 

One can then use two photon interference to find 
<pi j (26). If one inputs two photons into modes 
ii,%2 and detect in modes j\,j2, then we have 
\S) = S^otjO) and |D = «|0). If single, 



indistinguishable photons are used, then the proba- 
bility of post-selecting this output is 



Pi 



indist 



Per(A (S ' T: 

A (S,T).(S,T) .(S,T) A (S,T) 
Iv 22 IV 12 IV 21 

{ T il,ji T i2,32) + ( T h J2 r «2 



X cos 



J n ,32 



y '2 



'12,32 1 



(18) 



If the photons launched into the individual modes 
are distinguishable, then we get the incoherent sum 
of their individual statistics. We again take a matrix 
permanent, but as this is an incoherent process one 
finds 



P 



dist 



( T «l,il T «2 ,32) + ( T «l ,i2 T «2,jl) 



(19) 



One can then perform a Hong-Ou-Mandel interfer- 
ence experiment and find that the resulting dip vis- 
ibility is 



V 



P 



indist 



Pd, 



ist 



Pindist 

^ T h ,ji T h ,32 T i2 ,31 T i2 ,32 



(20) 



( T il,jl T i2,j2) 2 "I" i T ii,j2 T i2,ji) 2 
X COs{(j)i ujl + (j) i2 j 2 - 4> h j 2 - (j) i2jl ) (21) 

By measuring the dip visibility, and using the 
known one can find \4>a,ji + 4>i2,j2 - 4>n,j2 - 
<Pi2,ji\- Repeating this procedure for all accessi- 
ble two photon dips, and applying additional con- 
straints one can determine fa j (26). In theory, one 
can directly use Eqn. (18) to constrain the phases 
similarly to Eqn. (21). However, the former rela- 
tion is useful only if one post-selects on the entire 
N photon boson distribution which, in our case, re- 
quires number-resolving detectors. For our circuit 
topology, we encounter an overconstrained prob- 
lem as we measure more interference visibilities 
than unknown <p. Therefore, we run a least squares 
minimization to find the set of 4> that best fits our 
two photon interference data. This characterization 
process is efficient, as a general linear transforma- 
tion over M modes can be described by 0(M 2 ) 
parameters, requiring 0(M 2 ) measurements with 
this technique. 
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With Ajj experimentally determined over the 
accessible modes, one can predict the post-selected 
boson distribution for any input/output | S) , \T) by 
constructing A^ S ' T ) from A and taking the perma- 
nent according to Eqn. (13). Thus, we effectively 
use the one and two photon boson distributions to 
characterize our non-unitary operation over the ac- 
cessible modes. One can then take various matrix 
permanents of this non-unitary matrix to predict the 
boson distribution for any N. 

1.4 Sources of error 

The computational difficulty for a classical ma- 
chine to even verify the accuracy of QBSM output 
scales exponentially in N, making it vital to check 
the accuracy of these devices at low N where clas- 
sical computer validation remains easy. In this pa- 
per, we have benchmarked the accuracy of our QB- 
SMs by inferring the probability distribution from 
our data, labeled P cx p, and comparing it to the dis- 
tribution obtained from Equation (1), labeled P th . 
Throughout the text, we quantify the distance be- 
tween two probability distributions via the L\ dis- 
tance, d(pM p( 2 )) = \ £• |P; (1) - P ; (2) |. 

Our method of benchmarking QBSM accuracy 
will always yield a nonzero d due to the finite num- 
ber of samples collected. We perform a Monte 
Carlo simulation of a QBSM that perfectly sam- 
ples P th as a function of the number of counts col- 
lected, shown in Figure 4(a)-(b), to show the rate 
at which d asymptotically approaches zero as the 
number of samples collected increases. In the inset 
histograms, we show the range of outputs for this 
ideal QBSM for the actual number of experimental 
counts collected in the three and four photon cases, 
while the red dots indicate d(P th ,P exp ). These 
two probability distributions show close agreement 
in Figure 3, however a comparison of the histogram 
and red dots in Figure 4(a)-(b) indicates there is an 
additional source of error beyond the finite number 
of samples. 

Due to experimental limitations, occasionally we 
sample distributions other than P th . We model 
the effect of two such imperfections, photon dis- 
tinguishability and higher order terms from our 
PDC sources, and form a new distribution p mod 
that accounts for these effects. We find that the 



new <i(P exp ,P m ), indicated with the green dot 
in 4(a)-(b), comes well within the output variance 
of an ideal machine. This indicates we have cor- 
rectly diagnosed and modeled the principal sources 
of experimental error, which will be important in 
guiding designs of future, larger N, QBSMs. 

1.4.1 Effect of using heralded single photon 
sources 

The parametric downconversion sources we use to 
generate our photons actually generate a two-mode 
squeezed state that is given by 

oo 

I^pdc) = \A-A 2 ^A>n) (22) 

n=0 

where < A < 1 is the squeezing parameter whose 
magnitude is determined, in part, by the type of 
crystal and pump power used. For our experiment 
we wish to minimize higher order terms (ex. A 2 1 22) 
for the three photon QBSM) and so we deliberately 
lower our pump power as much as possible while 
maintaining a feasible count rate. For the three- 
photon experiment, A = VO.011 and for the four- 
photon experiment A = V0.023. However, even 
in this case we will sometimes inject more than A?~ 
photons into our circuit which, due to losses, could 
be observed as an A^-fold detection at the output. 

Due to our characterization method, we have no 
information about what such terms will be. In the 
operator sum representation, our characterization 
only determines the Eqq Kraus operator, which de- 
scribes the transformation of our input state when 
the environment, which includes all loss modes, 
starts and ends with zero photons. For example, 
if we instead inject five photons and lose two, then 
this process is described by the E20 operator, about 
which we have no information. 

To model the effect of using such squeezed 
sources, we start with a circuit with the same topol- 
ogy used in the experiment. Non-uniform losses 
throughout the circuit can then be modeled by 
adding beam splitters that link the depicted acces- 
sible modes shown in the figure to inaccessible loss 
modes, where the beam splitter reflectivity indi- 
cates the loss in that channel (33). Such 'loss beam- 
splitters' are added throughout the circuit. We can- 
not directly characterize these losses, however we 
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can estimate them numerically. The characterized 
linear transformation A is a function of these losses 
as well as the fabricated interferometers shown in 
Fig 2. We have performed a loss-independent char- 
acterization of the interferometers (20), and then 
input this data into a genetic algorithm to find the 
relative losses that best reproduces the character- 
ized A. 

With knowledge of these losses, we reproduce an 
actual unitary linear transformation U that extends 
over both the M = 6 modes as well as all loss 
modes and can be used to accurately predict the ef- 
fect of higher order PDC terms. We use Eqn. (1) 
to find the probability distributions when the input 
is |T) = 1011010), |022010) or |011020> which 
we label P 111 , p 221 and P 112 respectively. For ex- 
ample, the probability of obtaining |S) = 1 111000) 
given input |T) = 1 022010) is the summation of all 
terms where three photons appear in the first three 
accessible modes and two photons appear in any 
combination of loss modes. The higher order prob- 
ability distributions, p 221 and P 112 , are weighted 
by A 2 which is obtained via a conditional second 
order correlation measurement, g( 2 \0) (35). As 
A 2 is small, we only consider the first higher order 
terms from each source. 



|T) = |l) r |11000) where r labels a distinguish- 
able photon, then for a unitary transformation U 
the probability of obtaining an output |S) is 



P(S|T) 

+ 
+ 



U- 



U- 



U- 



(S,T) 

11 

(S,T) 
12 

(S,T) 
13 



\u. 



(S.T^S.T) 



22 



+ U. 



(S,T) [/ -(S,T)|2 



23 



32 



u. 



21 



33 



+ u. 



23 



31 



u. 



21 



32 



+ u. 



22 



31 



(24) 



where the terms in parentheses are permanents of 
2x2 matrices. We calculate these probability dis- 
tributions when one photon is distinguishable and 
weight them by |a 2 \/l — q 2 | 2 and |a 3 \/l — a 2 | 2 , 
the probability that one photon is distinguishable 
from the others for the two and three photon cases 
respectively. As a is large, we ignore the case when 
two photons are distinguishable. 



1.4.2 Photon Distinguishability 

Boson sampling assumes indistinguishable bosons, 
while experimental implementations will always 
have some distinguishability. Assuming pure in- 
puts we follow the notation of (28) to write our in- 
put state as 

N 

|V>in) = [J ( aA Li + V 7 ! - « 2 4,) |0) (23) 

i 

where N is again the number of photons, a is a dis- 
tinguishability parameter, and Al . f is the creation 
operator for photon i in mode £j. By analyzing the 
reduction in HOM dip visibility at a beamsplitter 
inside our circuit we estimate a = 0.974 in our 
experiment. 

If one photon is distinguishable from the others, 
then the new probability distribution is given by the 
permanents of iV — 1 matrices which are incoher- 
ently summed. For example, assume an input state 
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